Libraries

   if(any(grepl("package:dvutils", search()))) detach("package:dvutils", unload=TRUE)
   library(dvutils)
   library(stringr)
   library(ggplot2)
   library(dplyr)
   library(tidyr)
   library(tibble)
   library(readr)
   library(maps)
   library(glue)
   library(readxl)
   library(rlang)
   
   library(phyloseq)

   library(knitr)
   library(rmdformats)

   # library(pr2database)
   # data("pr2")

Markdown

Read the colors

class_colors.df <- read_excel("class_colors.xlsx", sheet = "classes")

class_colors <- structure(class_colors.df$color_name, .Names = class_colors.df$class)

Normalize and transform to long form

for (one_sample_type in sample_type_surface) {
    
    cat(one_sample_type)
    cat("\n============================")
    
    ps[[one_sample_type]] <- phyloseq_normalize_median(ps[[one_sample_type]])
    
    long[[one_sample_type]] <- phyloseq_transform_to_long(ps[[one_sample_type]])
    
    cat("\n")
    
}

18S filter surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 394 taxa and 50 samples ]
sample_data() Sample Data:       [ 50 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 394 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  23202
16S plastid surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 236 taxa and 40 samples ]
sample_data() Sample Data:       [ 40 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 236 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  28885
18S sort surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 114 taxa and 16 samples ]
sample_data() Sample Data:       [ 16 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 114 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  27275
18S filter 0.2 um surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 196 taxa and 17 samples ]
sample_data() Sample Data:       [ 17 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 196 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  14647
18S filter 3 um surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 229 taxa and 18 samples ]
sample_data() Sample Data:       [ 18 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 229 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  24783
18S filter 20 um surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 201 taxa and 15 samples ]
sample_data() Sample Data:       [ 15 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 201 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  29427
16S plastid 0.2 um surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 112 taxa and 11 samples ]
sample_data() Sample Data:       [ 11 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 112 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  28471
16S plastid 3 um surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 142 taxa and 16 samples ]
sample_data() Sample Data:       [ 16 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 142 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  31801
16S plastid 20 um surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 120 taxa and 13 samples ]
sample_data() Sample Data:       [ 13 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 120 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  25641
18S sort pico surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 65 taxa and 8 samples ]
sample_data() Sample Data:       [ 8 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 65 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  33258
18S sort nano surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 70 taxa and 8 samples ]
sample_data() Sample Data:       [ 8 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 70 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  21650

List of classes

classes <- data.frame(division = NULL, class = NULL)

for (one_sample_type in sample_type_surface[1:3]) {
    classes <- bind_rows(classes, select(long[[one_sample_type]], division, 
        class))
}

classes <- classes %>% distinct(division, class) %>% arrange(division, class)

# print(classes, n=50)
knitr::kable(classes)

Treemaps

treemap_class <- list()
treemap_genus <- list()

for (one_sample_type in sample_type_surface) {
    
    treemap_class[[one_sample_type]] <- phyloseq_long_treemap(long[[one_sample_type]], 
        division, class, str_c(one_sample_type, " - Class"), colors = class_colors)
    treemap_genus[[one_sample_type]] <- phyloseq_long_treemap(long[[one_sample_type]], 
        class, genus, str_c(one_sample_type, " - Genus"))
    
}

ASVs bargraphs

bargraph_asv <- list()
bargraph_species <- list()

for (one_sample_type in sample_type_surface) {
    
    bargraph_asv[[one_sample_type]] <- phyloseq_long_bargraph(long[[one_sample_type]], 
        text_scaling = 0.75, n_bars = 20, use_asv = TRUE, title = one_sample_type)
    bargraph_species[[one_sample_type]] <- phyloseq_long_bargraph(long[[one_sample_type]], 
        text_scaling = 0.75, n_bars = 20, use_asv = FALSE, title = one_sample_type)
    
}

## Barplots per sample

bargraph_sample <- list()

for (one_sample_type in sample_type_surface[4:11]) {
    
    bargraph_sample[[one_sample_type]] <- plot_bar(ps[[one_sample_type]], x = "sample_label", 
        fill = "class") + geom_bar(aes(color = class, fill = class), stat = "identity", 
        position = "stack") + ggtitle(str_c("Class level - ", one_sample_type)) + 
        theme(axis.text.y = element_text(size = 10)) + theme(axis.text.x = element_text(angle = 0, 
        hjust = 0.5)) + coord_flip() + scale_fill_viridis_d() + scale_color_viridis_d()
    print(bargraph_sample[[one_sample_type]])
}

Heatmaps

Filter : Dictyo / Chrysophyceae /Pelago /Bacili / Crypto / Pyrami / Mamiello

heatmap_class_abundant <- list()

for (one_sample_type in sample_type_surface[4:11]) {
    
    ps_heat <- tax_glom(ps[[one_sample_type]], taxrank = "class") %>% phyloseq_filter_abundant_taxa(fraction_min = 0.1)
    
    heatmap_class_abundant[[one_sample_type]] <- plot_heatmap(ps_heat, method = "NMDS", 
        distance = "bray", taxa.label = "class", taxa.order = "division", sample.label = "sample_label", 
        sample.order = "sample_label", low = "beige", high = "red", na.value = "beige", 
        trans = scales::log_trans(10), title = one_sample_type) + xlab("") + 
        ylab("")
    
    # plot(heatmap(otu_table(ps_heat)))
    
    print(heatmap_class_abundant[[one_sample_type]])
    
}

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 7 taxa and 17 samples ]
sample_data() Sample Data:       [ 17 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 7 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  13808

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 7 taxa and 18 samples ]
sample_data() Sample Data:       [ 18 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 7 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  23366

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 5 taxa and 15 samples ]
sample_data() Sample Data:       [ 15 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 5 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  27860

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 4 taxa and 11 samples ]
sample_data() Sample Data:       [ 11 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 4 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  27592

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 5 taxa and 16 samples ]
sample_data() Sample Data:       [ 16 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 5 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  31460

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 4 taxa and 13 samples ]
sample_data() Sample Data:       [ 13 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 4 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  24881

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 4 taxa and 8 samples ]
sample_data() Sample Data:       [ 8 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 4 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  32016

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 4 taxa and 8 samples ]
sample_data() Sample Data:       [ 8 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 4 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  21054

class_selected <- c("Dictyochophyceae", "Chrysophyceae", "Pelagophyceae", "Bacillariophyta", 
    "Cryptophyceae", "Pyramimonadales", "Mamiellophyceae")

heatmap_class_selected <- list()

for (one_sample_type in sample_type_surface[4:11]) {
    
    reads_max <- case_when(str_detect(one_sample_type, "18S filter") ~ 30000, 
        str_detect(one_sample_type, "16S plastid") ~ 30000, str_detect(one_sample_type, 
            "18S sort") ~ 25000)
    
    ps_heat <- tax_glom(ps[[one_sample_type]], taxrank = "class") %>% subset_taxa(class %in% 
        class_selected)
    
    # Try to order by division and class...  tax_table <-
    # data.frame(tax_table(ps_heat)@.Data) taxa_names(ps_heat) <-
    # str_c(tax_table$division, tax_table$class, sep=' - ')
    
    gg <- plot_heatmap(ps_heat, method = "NMDS", distance = "bray", taxa.label = "class", 
        taxa.order = "class", sample.label = "sample_label", sample.order = "sample_label", 
        low = "beige", high = "red", na.value = "gray95", trans = NULL, title = one_sample_type) + 
        xlab("") + ylab("") + theme(axis.text.x = element_text(angle = 45, hjust = 1, 
        vjust = 1)) + scale_fill_gradient(limits = c(0, reads_max), low = "beige", 
        high = "red", na.value = "gray95")
    
    if (str_detect(one_sample_type, "18S filter 0.2")) 
        gg <- gg + geom_vline(xintercept = c(6.5, 9.5, 12.5))
    
    if (one_sample_type == "18S filter 3 um surface") 
        gg <- gg + geom_vline(xintercept = c(8.5, 10.5, 13.5))
    
    if (one_sample_type == "18S filter 20 um surface") 
        gg <- gg + geom_vline(xintercept = c(8.5, 9.5))
    
    
    # plot(heatmap(otu_table(ps_heat))) , trans = scales::log_trans(10)
    
    print(gg)
    
    heatmap_class_selected[[one_sample_type]] <- gg
    
}

heatmap_genus_abundant <- list()

for (one_sample_type in sample_type_surface[4:11]) {
    
    ps_heat <- tax_glom(ps[[one_sample_type]], taxrank = "species") %>% phyloseq_filter_abundant_taxa(fraction_min = 0.1)
    
    heatmap_genus_abundant[[one_sample_type]] <- plot_heatmap(ps_heat, method = "NMDS", 
        distance = "bray", taxa.label = "species", taxa.order = "division", 
        sample.label = "sample_label", sample.order = "sample_label", low = "beige", 
        high = "red", na.value = "beige", title = one_sample_type)
    
    print(heatmap_genus_abundant[[one_sample_type]])
    
}

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 11 taxa and 17 samples ]
sample_data() Sample Data:       [ 17 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 11 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  11704

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 8 taxa and 18 samples ]
sample_data() Sample Data:       [ 18 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 8 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  20196

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 11 taxa and 15 samples ]
sample_data() Sample Data:       [ 15 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 11 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  26302

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 5 taxa and 11 samples ]
sample_data() Sample Data:       [ 11 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 5 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  27286

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 8 taxa and 16 samples ]
sample_data() Sample Data:       [ 16 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 8 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  30697

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 6 taxa and 13 samples ]
sample_data() Sample Data:       [ 13 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 6 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  24093

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 5 taxa and 8 samples ]
sample_data() Sample Data:       [ 8 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 5 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  28373

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 6 taxa and 8 samples ]
sample_data() Sample Data:       [ 8 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 6 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  19513

species_selected <- c("Thalassiosira_minima", "Fragilariopsis_cylindrus", "Minidiscus_sp.", 
    "Chaetoceros_neogracilis", "Porosira_glacialis", "Coretrhon_inerme", "Palmaria_palmata", 
    "Pseudo-nitzchia_seriata", "Micromonas_polaris", "Micromonas_clade_B3", 
    "Bathycoccus_prasinos", "Pyramimonas_gelidicola", "Geminigera_cryophila", 
    "Pelagophyceae_XXX_sp.", "Phaeocystis_antarctica")


heatmap_species_selected <- list()

for (one_sample_type in sample_type_surface[4:11]) {
    
    reads_max <- case_when(str_detect(one_sample_type, "18S filter") ~ 20000, 
        str_detect(one_sample_type, "16S plastid") ~ 5000, str_detect(one_sample_type, 
            "18S sort") ~ 15000)
    ps_heat <- tax_glom(ps[[one_sample_type]], taxrank = "species") %>% subset_taxa(species %in% 
        species_selected)
    
    # Try to order by division and species...  tax_table <-
    # data.frame(tax_table(ps_heat)@.Data) taxa_names(ps_heat) <-
    # str_c(tax_table$division, tax_table$species, sep=' - ')
    
    gg <- plot_heatmap(ps_heat, method = "NMDS", distance = "bray", taxa.label = "species", 
        taxa.order = "class", sample.label = "sample_label", sample.order = "sample_label", 
        low = "beige", high = "red", na.value = "gray95", trans = NULL, title = one_sample_type) + 
        xlab("") + ylab("") + theme(axis.text.x = element_text(angle = 45, hjust = 1, 
        vjust = 1)) + scale_fill_gradient(limits = c(0, reads_max), low = "beige", 
        high = "red", na.value = "gray95")
    
    if (str_detect(one_sample_type, "18S filter 0.2")) 
        gg <- gg + geom_vline(xintercept = c(6.5, 9.5, 12.5))
    
    if (one_sample_type == "18S filter 3 um surface") 
        gg <- gg + geom_vline(xintercept = c(8.5, 10.5, 13.5))
    
    if (one_sample_type == "18S filter 20 um surface") 
        gg <- gg + geom_vline(xintercept = c(8.5, 9.5))
    
    
    
    # plot(heatmap(otu_table(ps_heat))) , trans = scales::log_trans(10)
    
    print(gg)
    
    heatmap_species_selected[[one_sample_type]] <- gg
    
}

Apply to all samples

nmds_sample <- list()

for (one_sample_type in sample_type_surface[1:3]) {
    
    one_ps <- ps[[one_sample_type]] %>% tax_glom(taxrank = "species")
    
    nmds_sample[[one_sample_type]] <- phyloseq_nmds(one_ps, title = one_sample_type, 
        sample_color = season, sample_shape = fraction_name)
}

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1687566 
Run 1 stress 0.1687634 
... Procrustes: rmse 0.002985303  max resid 0.01690689 
Run 2 stress 0.1687197 
... New best solution
... Procrustes: rmse 0.002431247  max resid 0.01134347 
Run 3 stress 0.1733998 
Run 4 stress 0.1739068 
Run 5 stress 0.1766939 
Run 6 stress 0.1687159 
... New best solution
... Procrustes: rmse 0.001111716  max resid 0.005731252 
... Similar to previous best
Run 7 stress 0.1733984 
Run 8 stress 0.1739071 
Run 9 stress 0.1746154 
Run 10 stress 0.1835842 
Run 11 stress 0.168719 
... Procrustes: rmse 0.001095126  max resid 0.005701493 
... Similar to previous best
Run 12 stress 0.1835529 
Run 13 stress 0.168719 
... Procrustes: rmse 0.001110988  max resid 0.005739634 
... Similar to previous best
Run 14 stress 0.168719 
... Procrustes: rmse 0.001112726  max resid 0.005733226 
... Similar to previous best
Run 15 stress 0.1733994 
Run 16 stress 0.1745761 
Run 17 stress 0.1721531 
Run 18 stress 0.1734118 
Run 19 stress 0.1872138 
Run 20 stress 0.1687632 
... Procrustes: rmse 0.003837324  max resid 0.01748068 
*** Solution reached
[1] 3

Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance) 

global Multidimensional Scaling using monoMDS

Data:     wisconsin(sqrt(veganifyOTU(physeq))) 
Distance: bray 

Dimensions: 2 
Stress:     0.1687159 
Stress type 1, weak ties
Two convergent solutions found after 20 tries
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))' 

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1492484 
Run 1 stress 0.1492468 
... New best solution
... Procrustes: rmse 0.0002507616  max resid 0.001358323 
... Similar to previous best
Run 2 stress 0.1489429 
... New best solution
... Procrustes: rmse 0.01404688  max resid 0.06477591 
Run 3 stress 0.1732926 
Run 4 stress 0.1489447 
... Procrustes: rmse 0.0005970766  max resid 0.003273668 
... Similar to previous best
Run 5 stress 0.1526353 
Run 6 stress 0.1489494 
... Procrustes: rmse 0.001870841  max resid 0.01026137 
Run 7 stress 0.1489076 
... New best solution
... Procrustes: rmse 0.005462155  max resid 0.02114027 
Run 8 stress 0.1489094 
... Procrustes: rmse 0.0002241435  max resid 0.00122347 
... Similar to previous best
Run 9 stress 0.1931313 
Run 10 stress 0.1526371 
Run 11 stress 0.1669089 
Run 12 stress 0.2018991 
Run 13 stress 0.188055 
Run 14 stress 0.1526361 
Run 15 stress 0.1547389 
Run 16 stress 0.1760578 
Run 17 stress 0.1762028 
Run 18 stress 0.1965688 
Run 19 stress 0.1755688 
Run 20 stress 0.1870225 
*** Solution reached
[1] 3

Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance) 

global Multidimensional Scaling using monoMDS

Data:     wisconsin(sqrt(veganifyOTU(physeq))) 
Distance: bray 

Dimensions: 2 
Stress:     0.1489076 
Stress type 1, weak ties
Two convergent solutions found after 20 tries
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))' 

Square root transformation
Wisconsin double standardization
Run 0 stress 0.05090316 
Run 1 stress 0.05736207 
Run 2 stress 0.05090323 
... Procrustes: rmse 0.0001234125  max resid 0.0002263345 
... Similar to previous best
Run 3 stress 0.05185334 
Run 4 stress 0.05090317 
... Procrustes: rmse 3.166825e-05  max resid 5.733514e-05 
... Similar to previous best
Run 5 stress 0.05090316 
... Procrustes: rmse 3.131882e-05  max resid 5.132691e-05 
... Similar to previous best
Run 6 stress 0.05161249 
Run 7 stress 0.0520288 
Run 8 stress 0.05379784 
Run 9 stress 0.05723174 
Run 10 stress 0.05161249 
Run 11 stress 0.05379785 
Run 12 stress 0.05220006 
Run 13 stress 0.05940653 
Run 14 stress 0.05533946 
Run 15 stress 0.05344022 
Run 16 stress 0.05220002 
Run 17 stress 0.05161251 
Run 18 stress 0.05379788 
Run 19 stress 0.05161248 
Run 20 stress 0.0520288 
*** Solution reached
[1] 3

Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance) 

global Multidimensional Scaling using monoMDS

Data:     wisconsin(sqrt(veganifyOTU(physeq))) 
Distance: bray 

Dimensions: 2 
Stress:     0.05090316 
Stress type 1, weak ties
Two convergent solutions found after 20 tries
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))' 

Filter the data

sample_type_profile <- str_c(sample_type_all, " profile")

for (one_sample_type in sample_type_all) {
    
    one_sample_type_profile = str_c(one_sample_type, " profile")
    
    ps[[one_sample_type_profile]] <- ps[[one_sample_type]] %>% subset_samples(station_id == 
        "6") %>% subset_samples(date == "2015-01-16") %>% subset_samples(is.na(sample_concentration)) %>% 
        filter_taxa(function(x) sum(x) > 0, TRUE) %>% phyloseq_normalize_median()
    
    # Rename the samples
    
    if (!str_detect(one_sample_type_profile, "sort")) {
        sample_names(ps[[one_sample_type_profile]]) <- str_c(as.character(sample_data(ps[[one_sample_type_profile]])$depth), 
            "m", as.character(sample_data(ps[[one_sample_type_profile]])$fraction_min), 
            sep = "_")
    } else {
        sample_names(ps[[one_sample_type_profile]]) <- str_c(as.character(sample_data(ps[[one_sample_type_profile]])$depth), 
            "m", as.character(sample_data(ps[[one_sample_type_profile]])$fraction_name), 
            sep = "_")
    }
    
    
    print(glue("Phyloseq - {one_sample_type_profile}"))
    print(ps[[one_sample_type_profile]])
    cat("============================\n")
    
}

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 175 taxa and 15 samples ]
sample_data() Sample Data:       [ 15 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 175 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  43491Phyloseq - 18S filter profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 175 taxa and 15 samples ]
sample_data() Sample Data:       [ 15 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 175 taxa by 8 taxonomic ranks ]
============================

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 107 taxa and 15 samples ]
sample_data() Sample Data:       [ 15 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 107 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  41130Phyloseq - 16S plastid profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 107 taxa and 15 samples ]
sample_data() Sample Data:       [ 15 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 107 taxa by 8 taxonomic ranks ]
============================

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 96 taxa and 10 samples ]
sample_data() Sample Data:       [ 10 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 96 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  29305Phyloseq - 18S sort profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 96 taxa and 10 samples ]
sample_data() Sample Data:       [ 10 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 96 taxa by 8 taxonomic ranks ]
============================

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 92 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 92 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  30265Phyloseq - 18S filter 0.2 um profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 92 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 92 taxa by 8 taxonomic ranks ]
============================

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 86 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 86 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  31205Phyloseq - 18S filter 3 um profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 86 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 86 taxa by 8 taxonomic ranks ]
============================

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 121 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 121 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  58797Phyloseq - 18S filter 20 um profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 121 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 121 taxa by 8 taxonomic ranks ]
============================

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 46 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 46 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  40267Phyloseq - 16S plastid 0.2 um profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 46 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 46 taxa by 8 taxonomic ranks ]
============================

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 62 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 62 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  46172Phyloseq - 16S plastid 3 um profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 62 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 62 taxa by 8 taxonomic ranks ]
============================

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 74 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 74 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  47020Phyloseq - 16S plastid 20 um profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 74 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 74 taxa by 8 taxonomic ranks ]
============================

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 65 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 65 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  29412Phyloseq - 18S sort pico profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 65 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 65 taxa by 8 taxonomic ranks ]
============================

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 46 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 46 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  28502Phyloseq - 18S sort nano profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 46 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 46 taxa by 8 taxonomic ranks ]
============================

Barplots per depth

bargraph_sample <- list()

for (one_sample_type in sample_type_profile[4:11]) {
    
    bargraph_sample[[one_sample_type]] <- plot_bar(ps[[one_sample_type]], x = "depth", 
        fill = "class") + geom_bar(aes(color = class, fill = class), stat = "identity", 
        position = "stack") + ggtitle(str_c("Class level - ", one_sample_type)) + 
        theme(axis.text.y = element_text(size = 10)) + theme(axis.text.x = element_text(angle = 0, 
        hjust = 0.5)) + coord_flip() + # scale_fill_viridis_d() + scale_color_viridis_d() +
    scale_fill_manual(values = class_colors) + scale_color_manual(values = class_colors) + 
        scale_x_reverse() + xlab("Depth (m)") + ylab("Reads") + theme_bw()
    
    print(bargraph_sample[[one_sample_type]])
}

Sample by sample comparison

• Each dot corresponds to the relative contribution of the taxonomic level considered for a sample for which both 18S and 16S has been performed

long_species_sample <- list()
long_class_sample <- list()

for (one_sample_type in sample_type) {
    ps_one <- phyloseq_normalize_percent(ps[[one_sample_type]])
    
    long_one <- phyloseq_transform_to_long(ps_one) %>% group_by(date, station_id, 
        depth, fraction_name, division, class, order, family, genus, species) %>% 
        summarize(n_reads = sum(n_reads)) %>% ungroup()
    
    long_species_sample[[one_sample_type]] <- long_one
    
    long_one <- phyloseq_transform_to_long(ps_one) %>% group_by(date, station_id, 
        depth, fraction_name, division, class) %>% summarize(n_reads = sum(n_reads)) %>% 
        ungroup()
    
    long_class_sample[[one_sample_type]] <- long_one
}

long_compare_species_18S_16S <- full_join(rename(long_species_sample[["18S filter"]], 
    n_reads_18S_filter = n_reads), rename(long_species_sample[["16S plastid"]], 
    n_reads_16S_plastid = n_reads))

ggplot(long_compare_species_18S_16S, aes(x = n_reads_18S_filter, y = n_reads_16S_plastid)) + 
    geom_point(aes(color = class)) + scale_color_manual(values = class_colors) + 
    ggtitle("Species")

long_compare_class_18S_16S <- long_compare_species_18S_16S %>% group_by(date, 
    station_id, depth, fraction_name, division, class) %>% summarize(n_reads_18S_filter = sum(n_reads_18S_filter, 
    na.rm = TRUE), n_reads_16S_plastid = sum(n_reads_16S_plastid, na.rm = TRUE))

ggplot(long_compare_class_18S_16S, aes(x = n_reads_18S_filter, y = n_reads_16S_plastid)) + 
    geom_point(aes(color = class)) + scale_color_manual(values = class_colors) + 
    ggtitle("Classes")

Define upset function

• Needs only to have 0 and 1 in the matrix

plot_upset_long <- function(long_comparison, file_name) {
    
    long_upset <- long_comparison %>% mutate_if(is.numeric, funs(case_when(is.na(.) ~ 
        0, TRUE ~ 1))) %>% as.data.frame()
    
    fig_upset <- UpSetR::upset(long_upset, empty.intersections = "on", number.angles = 30, 
        point.size = 5, line.size = 2, mainbar.y.label = "Species intersections", 
        sets.x.label = "Species #", order.by = "freq", text.scale = c(2, 2, 
            2, 2, 2, 2), mb.ratio = c(0.7, 0.3))
    
    # text.scale = c(intersection size title, intersection size tick labels, set
    # size title, set size tick labels, set names, numbers above bars)
    
    print(fig_upset)
    png(file_name, width = 1000, height = 600, units = "px")
    print(fig_upset)
    dev.off()
}

Global comparison of species found with the 3 methods

• Compute table of number of samples for each species (rows) vs the three methods (columns)

• Only keep

  • 2015 samples because it is the only dataset for which we have the three types of samples
  • 0.2 and 3 um fractions (to be comparable with sorting)

• Remove species that contains _X

long_species <- list()
long_class <- list()

for (one_sample_type in sample_type) {
    long_one <- long_species_sample[[one_sample_type]] %>% filter(date >= as.Date("2015-01-01")) %>% 
        filter(!str_detect(species, "_X")) %>% filter(fraction_name != "> 20 um") %>% 
        group_by(division, class, order, family, genus, species) %>% summarize(n_reads = sum(n_reads), 
        n_samples = n()) %>% rename_at(vars(contains("n_")), funs(str_c(., str_replace_all(one_sample_type, 
        " ", "_"), sep = "_"))) %>% ungroup()
    
    long_species[[one_sample_type]] <- long_one
    
    long_one <- long_class_sample[[one_sample_type]] %>% filter(date >= as.Date("2015-01-01")) %>% 
        filter(fraction_name != "> 20 um") %>% group_by(division, class) %>% 
        summarize(n_reads = sum(n_reads), n_samples = n()) %>% rename_at(vars(contains("n_")), 
        funs(str_c(., str_replace_all(one_sample_type, " ", "_"), sep = "_"))) %>% 
        ungroup()
    
    long_class[[one_sample_type]] <- long_one
    
}

long_class_comparison <- long_class %>% purrr::reduce(full_join)
long_species_comparison <- long_species %>% purrr::reduce(full_join)

cat("Class not found in one type of sample\n")

Class not found in one type of sample

cat("18S filter\n")

18S filter

kable(filter(long_class_comparison, is.na(n_samples_18S_filter)))

cat("18S sort\n")

18S sort

kable(filter(long_class_comparison, is.na(n_samples_18S_sort)))

cat("16S filter\n")

16S filter

kable(filter(long_class_comparison, is.na(n_samples_16S_plastid)))

cat("Species  found in the three types of samples\n")

Species  found in the three types of samples

cat("18S filter\n")

18S filter

kable(filter(long_species_comparison, (!is.na(n_samples_18S_filter) & !is.na(n_samples_18S_sort) & 
    !is.na(n_samples_16S_plastid))))

cat("Species only found in one type of sample\n")

Species only found in one type of sample

cat("18S filter\n")

18S filter

kable(filter(long_species_comparison, (is.na(n_samples_18S_sort) & is.na(n_samples_16S_plastid))))

cat("18S sort\n")

18S sort

kable(filter(long_species_comparison, (is.na(n_samples_18S_filter) & is.na(n_samples_16S_plastid))))

cat("16S filter\n")

16S filter

kable(filter(long_species_comparison, (is.na(n_samples_18S_filter) & is.na(n_samples_18S_sort))))

openxlsx::write.xlsx(list(long_class_comparison, long_species_comparison), file = "../dada2/method_comparison.xlsx", 
    sheetName = c("Class", "Species"))

long_species_comparison_upset <- long_species_comparison %>% select(species, 
    contains("samples")) %>% rename_at(vars(contains("n_")), funs(str_replace(., 
    "n_samples_", "")))

plot_upset_long(long_species_comparison_upset, "fig/upset_method.png")

png 
  2 

Comparison of species found in the different fractions for 18S filter and 18S sorted

• Compute table of number of samples for each species (rows) vs the three fractions (columns)

• Only keep

  • 2015 samples because it is the only dataset for which we have the three types of samples
  • Genera that do not contain _X

long_species_fraction <- list()
long_class_fraction <- list()

for (one_sample_type in sample_type[c(1, 3)]) {
    long_one <- long_species_sample[[one_sample_type]] %>% filter(date >= as.Date("2015-01-01")) %>% 
        filter(!str_detect(species, "_X")) %>% group_by(fraction_name, division, 
        class, order, family, genus, species) %>% summarize(n_samples = n()) %>% 
        pivot_wider(names_from = fraction_name, values_from = n_samples) %>% 
        ungroup()
    
    long_species_fraction[[one_sample_type]] <- long_one
    
    long_one <- long_class_sample[[one_sample_type]] %>% filter(date >= as.Date("2015-01-01")) %>% 
        group_by(fraction_name, division, class) %>% summarize(n_samples = n()) %>% 
        pivot_wider(names_from = fraction_name, values_from = n_samples) %>% 
        ungroup()
    
    long_class_fraction[[one_sample_type]] <- long_one
    
}

long_class_fraction_comparison <- long_class_fraction %>% purrr::reduce(full_join)
long_species_fraction_comparison <- long_species_fraction %>% purrr::reduce(full_join)

kable(long_class_fraction_comparison)

kable(long_species_fraction_comparison)

plot_upset_long(long_species_fraction_comparison, "fig/upset_filter_sort.png")

png 
  2 

Comparison of species found in the different fractions for 18S filter

• Compute table of number of samples for each species (rows) vs the three fractions (columns)

• Only keep

  • Genera that do not contain _X

• We consider now ALL samples (will need to remove samples for which the > 20 um is missing)

long_species_fraction <- list()
long_class_fraction <- list()

for (one_sample_type in sample_type[c(1)]) {
    long_one <- long_species_sample[[one_sample_type]] %>% filter(!str_detect(species, 
        "_X")) %>% group_by(fraction_name, division, class, order, family, genus, 
        species) %>% summarize(n_samples = n()) %>% pivot_wider(names_from = fraction_name, 
        values_from = n_samples) %>% ungroup()
    
    long_species_fraction[[one_sample_type]] <- long_one
    
    long_one <- long_class_sample[[one_sample_type]] %>% group_by(fraction_name, 
        division, class) %>% summarize(n_samples = n()) %>% pivot_wider(names_from = fraction_name, 
        values_from = n_samples) %>% ungroup()
    
    long_class_fraction[[one_sample_type]] <- long_one
    
}

long_class_fraction_comparison <- long_class_fraction[[one_sample_type]]

long_species_fraction_comparison <- long_species_fraction[[one_sample_type]]

kable(long_class_fraction_comparison)

kable(long_species_fraction_comparison)

plot_upset_long(long_species_fraction_comparison, "fig/upset_filter.png")

png 
  2 

Comparison of species found in the different fractions for 18S sort

• Compute table of number of samples for each species (rows) vs the three fractions (columns)

• Only keep

  • Genera that do not contain _X

• We consider now ALL samples (will need to remove samples for which the > 20 um is missing)

long_species_fraction <- list()
long_class_fraction <- list()

for (one_sample_type in sample_type[c(3)]) {
    long_one <- long_species_sample[[one_sample_type]] %>% filter(!str_detect(species, 
        "_X")) %>% group_by(fraction_name, division, class, order, family, genus, 
        species) %>% summarize(n_samples = n()) %>% pivot_wider(names_from = fraction_name, 
        values_from = n_samples) %>% ungroup()
    
    long_species_fraction[[one_sample_type]] <- long_one
    
    long_one <- long_class_sample[[one_sample_type]] %>% group_by(fraction_name, 
        division, class) %>% summarize(n_samples = n()) %>% pivot_wider(names_from = fraction_name, 
        values_from = n_samples) %>% ungroup()
    
    long_class_fraction[[one_sample_type]] <- long_one
    
}

long_class_fraction_comparison <- long_class_fraction[[one_sample_type]]

long_species_fraction_comparison <- long_species_fraction[[one_sample_type]]

kable(long_class_fraction_comparison)

kable(long_species_fraction_comparison)